/* ======Trend Inflation and the Nature of Structural Breaks in the New Keynesian Phillips Curve ======

This program estimates the unobserved components model for the NKPC with stochastic trend inflation
with two-state Markov-switching. The variability of shocks to trend inflation are constrained to be constant
across regimes.

Written by: Pym Manopimoke

Model:

Pi_t = Pi_t^* + k_{S_t} Sum_{j=0}^{inf} beta^j E_{t-1}(x_{t+j}) + z_t
Pi_t^* = Pi_{t-1}^* + v_t, v_t ~iid.N(0,sig^2_v)
x_t = phi_1 x_{t-1} + phi_2 x_{t-2} + u_t, u_t|D_t~iid.N(0,sig^2_{u,D_t})
z_t = psi_{S_t} z_{t-1} + epsilon_t, epsilon_t|S_t ~iid.N(0,sig^2_{epsilon,S_t})

u_t and epsilon_t are correlated with correlation coefficient rho.

S_t = {0,1}
Prob[ S(t)=i | S(t-1)=j] = p_ij for i,j = 1,2.

sig^2_{u,D_t} = sig^2_{u,0}(1-D_t) + sig^2_{u,1}D_t
D_t = {0 if t<= 1984Q3, 1 otherwise}

========================================================================================================*/

new;
library optmum,PGRAPH;
pqgwin auto;
cls;
 
load data[250,11]=data.prn;  
data=data[2:250,.];
t=rows(data);
@Set the sample to 1952:1-2007:3@
startsmpl=21;
endsmpl =243; 

qtr= data[startsmpl:endsmpl,1];  
/*Column 2 for GDP chained-weighted inflation, 7 for CPI inflation, 8 for PCE inflation, 
9 for GDP Deflator*/  
inf=data[startsmpl:endsmpl,8];  @PCE inflation@
xx=data[startsmpl:endsmpl,5];  @CBO output gap@
t=rows(inf~xx);

@Define dummy variable series D_t to account for the Great Moderation@
GM = 130; @1984:2@
dmat = ones(GM,1)|zeros(t-GM,1); 
 

@================= Initialize Global Variables============@

tvp_no=4; 
start= 11;
t1=75;t2=125;
   
prmtr_in= 2 |
          -2 |
          1.0310126 |
          0.45915495 |
          0.29570007 |
          1.0090050 |
          2.2200207 |
          -0.079288579 |
          1.4791923 |
          -4.3570819 |
          0.91213739 |
          6.1668916 |
          0.018944040 |
          0.021589511 ;
  
{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);
PRM_FNL=TRANS(xout);             @ Estimated  coefficients, constrained@
xout;
hout0=hessp(&lik_fcn,xout);
hout0;
hout=inv(hout0);
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@
 
{B_TT,PR_TT,PR_TL} =FILTER(XOUT);

sm2=SMOOTH2(pr_tt,pr_tl,prm_fnl[1:2,1]); @Pr[S_t|Y_T], Smoothed probabilities@

"== FINAL OUTPUT ===============================";
"Likelihood value is ";; -fout;
"Estimated parameters (s.e.) tstat are:"; 
prm_fnl~sd_fnl~prm_fnl./sd_fnl;
"===============================================";

@Produces figure A1 in the unpublished appendix@
_pdate=0;
_pltype = {6 3};
_plegctl={4 5 5 4.8};
_plegstr="Regime1\000Regime2";
 title("Smoothed probabilities");
 title("");
 XY(qtr[start:t],sm2[.,1]~sm2[.,3]);

_pltype = {6 3};
_plegctl={2 4 6 4.8};
_plegstr="Trend Inflation\000Actual Inflation";
title("Actual inflation and trend inflation");
xy(qtr[start:t],b_tt[.,1]~inf[start:t]);

output file=UC1.out reset;
output off;
end;

@ end OF MAIN PROGRAM @
@========================================================================@
 

PROC LIK_FCN(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,ppr,qpr,prob_1,prob_0,
xfor,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2, 
B_TL00,B_TL01,B_TL10,B_TL11,P_TL00,P_TL01,P_TL10,P_TL11,F_CAST00,F_CAST01,F_CAST10,F_CAST11,SS00,SS01,SS10,SS11,
B_LL0,B_LL1,P_LL0,P_LL1,pr_ll0_,pr_ll1_,B_TT00,B_TT01,B_TT10,B_TT11,P_TT00,P_TT01,P_TT10,P_TT11,
PR_VL00,PR_VL01,PR_VL10,PR_VL11,PR_VAL,PRO_00,PRO_01,PRO_10,PRO_11;  

EXTERNAL PROC TRANS, V_PROB;
 
           PRMTR=TRANS(PRMTR1);

            ppr=PRMTR[1,1];    @Pr[St=1/St-1=1]@
            qpr=PRMTR[2,1];    @Pr[St=0/St-1=0]@
            prob_1=(1-qpr)/(2-ppr-QPR); @Pr[St-1=1/Yt-1], Steady State Prob.@
            prob_0=1-PROB_1;  @Pr[St-1=1/Yt-1], Steady State Prob.@
         
            siguGM0 = prmtr[3,1];
            siguGM1 = prmtr[4,1];
            sigv= prmtr[5,1]; 
            sigep = prmtr[6:7,1];  
            rho = prmtr[8,1];
            phi1 = prmtr[9,1];
            phi2 = prmtr[10,1];
            psi=prmtr[11:12,1];    
            k=prmtr[13:14,1];

            A=1~0; @selection matrix@
            H = 1~1~0~0|0~0~1~0;
                         
            B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
            P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
            B_LL1=inf[t1]|0|xx[1]|0; 
            P_LL1 = eye(4)*10000;       
                
            subF = phi1~phi2| 1~0; 
            subQ= siguGM0^2~0|0~0;
            vecp = inv(eye(4)-subF.*.subf)*vec(subq);
            vv = reshape(vecp',2,2);
            p_ll0[3:4,3:4] = vv;
            p_ll1[3:4,3:4] = vv;
                  
            p_ll0[2,2] = sigep[1,1]^2/(1-psi[1]^2);
            p_ll1[2,2] = sigep[2,1]^2/(1-psi[2]^2);
 

        /*  START ITERATION */
        /*===================================================================== */
        
           LIKV=0.0;
           J_ITER = 1;
           do until J_ITER>T;

                yt=inf[j_iter]|xx[j_iter];
                dmt = dmat[j_iter];
                sigu=siguGM0*dmt + siguGM1*(1-dmt);
                
                Q0=    sigv^2~0~0~0|
                       0~sigep[1]^2~rho*sigu*sigep[1,1]~0|
                       0~rho*sigu*sigep[1,1]~sigu^2~0|
                       0~0~0~0;
                Q1=    sigv^2~0~0~0|
                       0~sigep[2]^2~rho*sigu*sigep[2,1]~0|
                       0~rho*sigu*sigep[2,1]~sigu^2~0|
                       0~0~0~0;
             
                F0 =   1~0~0~0 |
                       0~psi[1]~0~0|
                       0~0~phi1~phi2|
                       0~0~1~0;
                F1 =   1~0~0~0 |
                       0~psi[2]~0~0|
                       0~0~phi1~phi2|
                       0~0~1~0;
                    
                FF= phi1~phi2|1~0;
                xfor = k[1]*A*(inv(eye(2)-FF));  @infinite sum term in cyclical component@        
                ome0 =0~0~xfor| 0~0~0~0;
                xfor = k[2] *A*(inv(eye(2)-FF)); 
                ome1 =0~0~xfor| 0~0~0~0;
                    
                /* KALMAN FILTER */
                @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
                B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
                B_TL01 = F1 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
                B_TL10 = F0 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
                B_TL11 = F1 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@
                
                @ P_tl(i,j),  S(t-1)=i, S(t)=j@
                P_TL00 = F0 * P_LL0 * F0' + Q0;
                P_TL01 = F1 * P_LL0 * F1' + Q1;
                P_TL10 = F0 * P_LL1 * F0' + Q0;
                P_TL11 = F1 * P_LL1 * F1' + Q1;
            
                @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
                F_CAST00= yt- H * B_TL00 - ome0* B_TL00;
                F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
                F_CAST10= yt- H * B_TL10 - ome0* B_TL10;
                F_CAST11= yt- H * B_TL11 - ome1* B_TL11;
                    
                @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
                SS00= H * P_TL00 * H' ; 
                SS01= H*  P_TL01 * H' ;
                SS10= H * P_TL10 * H' ;
                SS11= H * P_TL11 * H' ;
              
                /*Updating*/
                B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
                B_TT01 = B_TL01 + (P_TL01 * H' *inv(SS01)) * F_CAST01;
                B_TT10 = B_TL10 + (P_TL10 * H' *inv(SS10)) * F_CAST10;
                B_TT11 = B_TL11 + (P_TL11 * H' *inv(SS11)) * F_CAST11;
             
                P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
                P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
                P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
                P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
               
                /* HAMILTON FILTER */   
                PR_VL00=V_PROB(F_CAST00,SS00)*qpr *PROB_0;    @Pr[St,Yt/Yt-1]@
                PR_VL01=V_PROB(F_CAST01,SS01)*(1-qpr)*PROB_0;
                PR_VL10=V_PROB(F_CAST10,SS10)*(1-ppr)*PROB_1;
                PR_VL11=V_PROB(F_CAST11,SS11)*ppr *PROB_1;
             
                PR_VAL=PR_VL00+PR_VL01+PR_VL10+PR_VL11;
    
                @Pr[St,St-1/Yt]@
                PRO_00=PR_VL00/PR_VAL;       
                PRO_01=PR_VL01/PR_VAL;
                PRO_10=PR_VL10/PR_VAL;
                PRO_11=PR_VL11/PR_VAL;
           
                PROB_0=PRO_00+PRO_10;       @Pr[St=0/Yt]@
                PROB_1=PRO_01+PRO_11;       @Pr[St=1/Yt]@

                /* COLLAPSING TERMS */
   
                B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10  )/PROB_0;
                B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11  )/PROB_1;
                        
                P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                         PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)'))/PROB_0;
                P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                         PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)'))/PROB_1;
                 
                if J_ITER < START; goto skip; endif;
 
                LIKV = LIKV -LN(PR_VAL);
                skip:
            J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
  
RETP(LIKV);
ENDP;



@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc TRANS(prmtr1); @ constraining values of reg. coeff.@
local prm_new,d0,d1;
prm_new=prmtr1; 

prm_new[1:2]=EXP(-prmtr1[1:2]) ./(1+exp(-1*prmtr1[1:2]));  @constraining the probabilities@ 
prm_new[3:7,1] = sqrt(prmtr1[3:7,1]^2);
prm_new[8] =prmtr1[8,1]./(1+abs(prmtr1[8,1]));
d0 =prmtr1[9,1]./(1+abs(prmtr1[9,1]));
d1 =(1-abs(d0))*prmtr1[10,1]./(1+abs(prmtr1[10,1]))+ abs(d0)-d0^2;
prm_new[9,1]=2*d0;
prm_new[10,1]=-1*(d0^2 + d1) ;
prm_new[11:12,1] =prmtr1[11:12,1]./(1+abs(prmtr1[11:12,1]));
 
retp(prm_new);
endp;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
PROC V_PROB(EV, HE);      @ CALCULATES    Pr[Yt/St,Yt-1]@
local VAL;

VAL=(1/SQRT(2*PI*det(HE)))*EXP(-0.5*EV'inv(HE)*EV);

RETP(VAL);
ENDP;
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@


proc (3) = FILTER(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,ppr,qpr,prob_1,prob_0,b_m, pr_ttm, pr_tlm,
xfor,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2,pr_0_l,pr_1_l,
B_TL00,B_TL01,B_TL10,B_TL11,P_TL00,P_TL01,P_TL10,P_TL11,F_CAST00,F_CAST01,F_CAST10,F_CAST11,SS00,SS01,SS10,SS11,
B_LL0,B_LL1,P_LL0,P_LL1,pr_ll0_,pr_ll1_,B_TT00,B_TT01,B_TT10,B_TT11,P_TT00,P_TT01,P_TT10,P_TT11,
PR_VL00,PR_VL01,PR_VL10,PR_VL11,PR_VAL,PRO_00,PRO_01,PRO_10,PRO_11;  

EXTERNAL PROC TRANS, V_PROB;
 
            PRMTR=TRANS(PRMTR1);
 
            ppr=PRMTR[1,1];    @Pr[St=1/St-1=1]@
            qpr=PRMTR[2,1];    @Pr[St=0/St-1=0]@
            prob_1=(1-qpr)/(2-ppr-QPR); @Pr[St-1=1/Yt-1], Steady State Prob.@
            prob_0=1-PROB_1;  @Pr[St-1=1/Yt-1], Steady State Prob.@
         
            siguGM0 = prmtr[3,1];
            siguGM1 = prmtr[4,1];
            sigv= prmtr[5,1]; 
            sigep = prmtr[6:7,1];  
            rho = prmtr[8,1];
            phi1 = prmtr[9,1];
            phi2 = prmtr[10,1];
            psi=prmtr[11:12,1];    
            k=prmtr[13:14,1];

            A=1~0; @selection matrix@
            H = 1~1~0~0|0~0~1~0;
                         
            B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
            P_LL0 = eye(4)*10000;          @uncertainty associated with initial guess@
            B_LL1=inf[t1]|0|xx[1]|0; 
            P_LL1 = eye(4)*10000;       
                
            subF = phi1~phi2| 1~0; 
            subQ= siguGM0^2~0|0~0;
            vecp = inv(eye(4)-subF.*.subf)*vec(subq);
            vv = reshape(vecp',2,2);
            p_ll0[3:4,3:4] = vv;
            p_ll1[3:4,3:4] = vv;
                  
            p_ll0[2,2] = sigep[1,1]^2/(1-psi[1]^2);
            p_ll1[2,2] = sigep[2,1]^2/(1-psi[2]^2);
  
            B_M=ZEROS(T,tvp_no);  @will store (1,1)element of B_TT@
            PR_TTM=ZEROS(T,2);  @will store Pr[S_t=0|Y_t]@
            PR_TLM=ZEROS(T,2);  @will store Pr[S_t=0|Y_{t-1}]@
 
        /*  START ITERATION */
 
            LIKV=0.0;
            J_ITER = 1;
            do until J_ITER>T;
               
                yt=inf[j_iter]|xx[j_iter];

                pr_0_L=qpr*prob_0+(1-ppr)*prob_1;  @PR[St=0/Yt-1]@
                pr_1_L=(1-QPR)*prob_0+ppr*prob_1;  @PR[St=1/Yt-1]@ 
                PR_TLM[J_ITER,.]= pr_0_L~pr_1_L;
         
                dmt = dmat[j_iter];
                sigu=siguGM0*dmt + siguGM1*(1-dmt);
                
                Q0=    sigv^2~0~0~0|
                       0~sigep[1]^2~rho*sigu*sigep[1,1]~0|
                       0~rho*sigu*sigep[1,1]~sigu^2~0|
                       0~0~0~0;
                Q1=    sigv^2~0~0~0|
                       0~sigep[2]^2~rho*sigu*sigep[2,1]~0|
                       0~rho*sigu*sigep[2,1]~sigu^2~0|
                       0~0~0~0;
             
                F0 =   1~0~0~0 |
                       0~psi[1]~0~0|
                       0~0~phi1~phi2|
                       0~0~1~0;
                F1 =   1~0~0~0 |
                       0~psi[2]~0~0|
                       0~0~phi1~phi2|
                       0~0~1~0;
                    
                FF= phi1~phi2|1~0;
                xfor = k[1]*A*(inv(eye(2)-FF));  @infinite sum term in cyclical component@        
                ome0 =0~0~xfor| 0~0~0~0;
                xfor = k[2] *A*(inv(eye(2)-FF)); 
                ome1 =0~0~xfor| 0~0~0~0;

                /* KALMAN FILTER */
                @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
                B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
                B_TL01 = F1 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
                B_TL10 = F0 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
                B_TL11 = F1 * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@
                
                @ P_tl(i,j),  S(t-1)=i, S(t)=j@
                P_TL00 = F0 * P_LL0 * F0' + Q0;
                P_TL01 = F1 * P_LL0 * F1' + Q1;
                P_TL10 = F0 * P_LL1 * F0' + Q0;
                P_TL11 = F1 * P_LL1 * F1' + Q1;
            
                @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
                F_CAST00= yt- H * B_TL00 - ome0* B_TL00;
                F_CAST01= yt- H * B_TL01 - ome1* B_TL01;
                F_CAST10= yt- H * B_TL10 - ome0* B_TL10;
                F_CAST11= yt- H * B_TL11 - ome1* B_TL11;
                    
                @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
                SS00= H * P_TL00 * H' ; 
                SS01= H*  P_TL01 * H' ;
                SS10= H * P_TL10 * H' ;
                SS11= H * P_TL11 * H' ;
              
                /*Updating*/
                B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
                B_TT01 = B_TL01 + (P_TL01 * H' *inv(SS01)) * F_CAST01;
                B_TT10 = B_TL10 + (P_TL10 * H' *inv(SS10)) * F_CAST10;
                B_TT11 = B_TL11 + (P_TL11 * H' *inv(SS11)) * F_CAST11;
             
                P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
                P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
                P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
                P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
            
                /* HAMILTON FILTER */   
                PR_VL00=V_PROB(F_CAST00,SS00)*qpr *PROB_0;    @Pr[St,Yt/Yt-1]@
                PR_VL01=V_PROB(F_CAST01,SS01)*(1-qpr)*PROB_0;
                PR_VL10=V_PROB(F_CAST10,SS10)*(1-ppr)*PROB_1;
                PR_VL11=V_PROB(F_CAST11,SS11)*ppr *PROB_1;
             
                PR_VAL=PR_VL00+PR_VL01+PR_VL10+PR_VL11;
    
                @Pr[St,St-1/Yt]@
                PRO_00=PR_VL00/PR_VAL;       
                PRO_01=PR_VL01/PR_VAL;
                PRO_10=PR_VL10/PR_VAL;
                PRO_11=PR_VL11/PR_VAL;
           
                PROB_0=PRO_00+PRO_10;       @Pr[St=0/Yt]@
                PROB_1=PRO_01+PRO_11;       @Pr[St=1/Yt]@

                /* COLLAPSING TERMS */
   
                B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10  )/PROB_0;
                B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11  )/PROB_1;
                        
                P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                         PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)'))/PROB_0;
                P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                         PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)'))/PROB_1;

                @Estimates of the unobserved components@
                B_TT=B_TT00*PRO_00+B_TT01*PRO_01+B_TT10*PRO_10+B_TT11*PRO_11;
                B_M[J_ITER,.]=B_TT';
                
                PR_TTM[J_ITER,.]=PROB_0~prob_1;
  
                if J_ITER < START; goto skip; endif;
                LIKV = LIKV -LN(PR_VAL);

                skip:
            J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
/*===================*/

retp(B_M[START:T,.],PR_TTM[START:T,.],PR_TLM[START:T,.]);
ENDP;


@==========================================================@

@Kim's smoothing@
proc SMOOTH2(pr_tt0,pr_tl0,prmtr0);

local ppr, qpr, pr_sm0,pr_sm1,j_iter,pr_sm00, pr_sm01, pr_sm10, pr_sm11,TT,prmtr,pr_tt1, pr_tl1  ;

     TT=ROWS(PR_TT0);
     prmtr=prmtr0;
 
     ppr = prmtr[1,1];
     qpr = prmtr[2,1];
     
     pr_tt1=1-pr_tt0;
     pr_tl1=1-pr_tl0;
 
     pr_sm0=pr_tt0;      
     pr_sm1=pr_tt1;
  
     j_iter=TT-1;
     do until j_iter < 1;
  
        pr_sm00=pr_sm0[j_iter+1,1]*qpr*pr_tt0[j_iter,1]/pr_tl0[j_iter+1,1];
        pr_sm01=pr_sm1[j_iter+1,1]*(1-qpr)*pr_tt0[j_iter,1]/pr_tl1[j_iter+1,1];
        pr_sm10=pr_sm0[j_iter+1,1]*(1-ppr)*pr_tt1[j_iter,1]/pr_tl0[j_iter+1,1];
        pr_sm11=pr_sm1[j_iter+1,1]*ppr*pr_tt1[j_iter,1]/pr_tl1[j_iter+1,1];
        
        pr_sm0[j_iter,1]=pr_sm00+pr_sm01;
        pr_sm1[j_iter,1]=pr_sm10+pr_sm11;
 
     j_iter=j_iter -1;
     endo;
 
retp(pr_sm0~pr_sm1);  
endp;

/*=================================================================================================================*/
 
 